# Author: Darren Colby
# Date: 10/6/2021
# Purpose: To examine the fragmentation of cartels and militias after
# El Chapo's capture

# Initial settings --------------------------------------------------------

library(blockmodeling)
library(tidyverse)
library(ggraph)
library(igraph)
library(ggthemes)
library(patchwork)
library(texreg)
library(RSiena)
library(extrafont)

# Load helper functions
source("scripts/helper_functions.R")

# Load fonts for making graphs
loadfonts(device = "win")

# Load the data -----------------------------------------------------------

drugnet <- read_csv("output/drugnet.csv")

edge_df <- drugnet %>% 
   select(sideA, sideB, weight)

combined_graph <- graph_from_data_frame(edge_df)

combined_matrix <- as.matrix(as_adjacency_matrix(combined_graph))

main_graph <- edge_df %>% 
   distinct(sideA, sideB) %>% 
   graph_from_data_frame()

# Create roles for each actor ---------------------------------------------

# Loop through each number of equivalence groups and find the mean error
block_errors <- map(seq(2, 10, 1), function(k){
   mean(optRandomParC(M = combined_matrix, k = k, rep = 10, approaches = "ss",
                 blocks = "com", seed = 28)$err)}) %>% 
   unlist()

# Plot the structural equivalence solutions -------------------------------

# Create a dataframe from the blockmodels
bind_cols(block_errors, seq(2, 10, 1)) %>% 
   rename("mean_error" = ...1, "groups" = ...2) %>% 
   
   # Plot it
   ggplot(aes(x = groups,
           y = mean_error)) + 
   geom_vline(aes(xintercept = 7),
              linetype = "dashed",
              color = "gray") + 
   geom_hline(aes(yintercept = 56698.39),
              linetype = "dashed",
              color = "gray") + 
   geom_line() + 
   geom_point() + 
   scale_x_continuous(breaks = seq(2, 10, 1),
                      labels = c("2", "3", "4", "5", "6", "Best solution", 
                                 "8", "9", "10")) + 
   labs(x = "Groups",
        y = "Mean error",
        caption = "Note: 10 iterations were used for each number of partitions") + 
   theme_few()

ggsave("figures/figureA1.tiff", width = 6.5, height = 4, dpi = 1000)

# Add structural equivalence group attribute ------------------------------

# Create a 7 block solution
blocks <- optRandomParC(M = combined_matrix, 
                        k = 7,  
                        rep = 100,  
                        approaches = "ss", 
                        blocks = "com",  
                        seed = 42)

# Assign actors to blocks
V(main_graph)$role <- blocks$best$best1$clu

# Convert roles to names
V(main_graph)$role <- case_when(
   V(main_graph)$role == 1 ~ "Gulf Cartel",
   V(main_graph)$role == 2 ~ "Rising Challengers",
   V(main_graph)$role == 3 ~ "Sinaloa Cartel",
   V(main_graph)$role == 4 ~ "Small Cartels and Militias",
   V(main_graph)$role == 5 ~ "Los Zetas",
   V(main_graph)$role == 6 ~ "White Dwarfs",
   V(main_graph)$role == 7 ~ "Cartel Jalisco Nueva Generacion",
   TRUE ~ as.character(V(main_graph)$role)
)

# Add attribute for militia status ----------------------------------------

militias <- c("malinaltepec communal militia (mexico)", "chamula militia", 
              "san pablo cuatro venados communal militia (mexico)", 
              "dzan communal militia (mexico)", 
              "el pinar communal militia (mexico)", 
              "loma de la cruz communal militia (mexico)", 
              "molinos los arcos communal militia (mexico)", 
              "el carrizal de bravo communal militia (mexico)", 
              "huahua communal militia (mexico)", 
              "tlacotepec communal militia (mexico)", 
              "maya indigenous militia (mexico)", 
              "mixe indigenous militia (mexico)", 
              "san agustin oapan communal militia (mexico)", 
              "santa isabel de la reforma communal militia (mexico)", 
              "aldama indigenous militia (mexico)", 
              "san mateo yucutindoo communal militia (mexico)", 
              "apetlanca communal militia (mexico)", 
              "rio santiago communal militia (mexico)",
              "yautepec communal militia (mexico)", 
              "upoeg self defense group", 
              "faction of front for security and development vigilante group", 
              "totolapan self-defense force", "autodefensas unidas de michoacan", 
              "alacatlatzala communal militia (mexico)", 
              "chocomanatlan communal militia (mexico)", 
              "coahuayutla de guerrero communal militia (mexico)", 
              "cotija de la paz communal militia (mexico)", 
              "cuilapam de guerrero communal militia (mexico)", 
              "el cipresal communal militia (mexico)", 
              "el nith communal militia (mexico)", 
              "ezln: zapatista army of national liberation", 
              "leonardo bravo communal militia (mexico)", 
              "los reyes communal militia (mexico)", 
              "san cristobal de las casas communal militia (mexico)", 
              "san miguel tecuiciapan communal militia (mexico)", 
              "santa maria nativitas coatlan communal militia (mexico)", 
              "santa martha indigenous militia (mexico)", 
              "santiago amoltepec communal militia (mexico)", 
              "santiago miltepec communal militia (mexico)", 
              "santo reyes zochiquilazala communal militia (mexico)", 
              "tangancicuaro communal militia (mexico)", 
              "tenquizolco communal militia (mexico)", 
              "tepalcatepec communal militia (mexico)", 
              "tinguindin communal militia (mexico)", 
              "tocumbo communal militia (mexico)", 
              "villa morelos communal militia (mexico)", 
              "villa victoria communal militia (mexico)", 
              "xochiltepec communal militia (mexico)",
              "rival faction of front for security and development vigilante group")

# Assign militia status to nodes
V(main_graph)$militia <- ifelse(V(combined_graph)$name %in% militias, 1, 0)

# Assign aggressiveness attribute -----------------------------------------

# Convert the graph to a dataframe
temp_vector <- as_data_frame(main_graph) %>% 
   as_tibble() %>% 
   select(from) %>% 
   
   # Merge the vertices with the the aggression data for sending actors
   left_join(drugnet %>% 
                select(sideA, agression) %>% 
                distinct(), 
             by = c("from" = "sideA")) %>% 
   distinct(from, agression) %>% 
   
   # Add the sideB actors to the dataframe. Using a full join will only add rows
   # for actors not already in the dataframe
   full_join(drugnet %>% 
                select(sideB, agression) %>% 
                distinct(),
             by = c("from" = "sideB")) %>% 
   select(-agression.y) %>% 
   rename("agression" = agression.x) %>% 
   distinct(from, agression) %>% 
   select(agression) %>% 
   
   # Codes actors that never directed an attack as having an aggression of 0
   mutate(agression = ifelse(is.na(agression), 0, agression)) %>% 
   ungroup() %>% 
   pull()

# Assign aggression scores to vertices
V(main_graph)$aggression <- temp_vector 

# Code subfaction status --------------------------------------------------

subfaction <- vector(length = 151, mode = "numeric")

# The indices of groups that are subfactions
subfaction[c(8, 72, 74, 75, 77:79, 113, 134, 136:138, 142:146, 148, 149)] <- 1
V(main_graph)$subfaction <- subfaction

# Create graphs for each time period --------------------------------------

# Create a vector of periods in which each tie existed
period <- drugnet %>% 
   mutate(period = ifelse(year < 2017, 1, 2)) %>% 
   group_by(period, sideA, sideB) %>% 
   ungroup() %>% 
   distinct(sideA, sideB, period) %>% 
   group_by(sideA, sideB) %>% 
   summarise(period = max(period)) %>% 
   ungroup() %>% 
   pull(period)

# Add a year attribute to the combined graph
E(main_graph)$period <- period

# Split into two graphs
pre_arrest_graph <- subgraph.edges(main_graph, 
                                   E(main_graph)[E(main_graph)$period < 2])

post_arrest_graph <- subgraph.edges(main_graph, 
                                    E(main_graph)[E(main_graph)$period > 1])

# Plot the graphs before and after El Chapo's arrest ----------------------

pre_arrest_plot <- plot_network(pre_arrest_graph, "Pre-arrest") + 
   theme(legend.position = "none")
post_arrest_plot <- plot_network(post_arrest_graph, "Post-arrest")

pre_arrest_plot / post_arrest_plot + plot_layout(guides = "collect")


ggsave("figures/figure1.tiff", width = 10, height = 8.5, dpi = 600)

# Plot degree distribution ------------------------------------------------

p3 <- plot_degree_distribution(pre_arrest_graph, "in", "Before", "In-degree", 
                               c(0, 25))
p4 <- plot_degree_distribution(pre_arrest_graph, "out", NULL, "Out-degree", 
                               c(0, 20))
p5 <- plot_degree_distribution(post_arrest_graph, "in", "After", "In-degree", 
                               c(0, 105))
p6 <- plot_degree_distribution(post_arrest_graph, "out", NULL, "Out-degree", 
                               c(0, 85))

p3 + p5 +  p4 + p6

ggsave("figures/figure2.tiff", width = 6.2, height = 3.7, dpi = 600)

# Prepare the data for SAOM estimation ------------------------------------

# Make a graph with only edges in the first period
t0_network <- main_graph %>% 
   delete_edges(as.vector(E(post_arrest_graph)))

# Make a graph with only edges in the second period
t1_network <- subgraph.edges(main_graph, as.vector(E(post_arrest_graph)), 
                             delete.vertices = FALSE)

# Make a list for joiner nodes
joiner_nodes <- as_tibble(as.vector(V(t0_network)$name)) %>% 
                   mutate(node = row_number()) %>% 
                   ungroup() %>% 
                   relocate(node) %>% 
                   filter(!(value %in% as.vector(V(pre_arrest_graph)$name))) %>% 
                   pull(node)

# Make a list for all nodes and code them as being present in both periods
joiners <- rep(list(c(1, 2)), 151)

# Update periods for nodes that were only present in the second period
for (i in seq_along(joiners)){
   
   # Check to see if a node was a joiner in the second period
   if (i %in% joiner_nodes){
      
      # Recode the list to show joiners only in period 2
      joiners[[i]][1] <- 2
      
   }
   
}

# Creates an RSiena object for joiners
joiners <- sienaCompositionChange(joiners)

# Create a dependent variable
dv <- sienaDependent(array(c(as.matrix(as_adjacency_matrix(t0_network)), 
                             as.matrix(as_adjacency_matrix(main_graph))), 
                           dim = c(151, 151, 2)))

# Create constant covariates objects
role_var <- coCovar(as.vector(blocks$best$best1$clu))
militia_var <- coCovar(as.vector(V(main_graph)$militia))
aggression_var <- coCovar(as.vector(V(main_graph)$aggression))
subfaction_var <- coCovar(as.vector(V(main_graph)$subfaction))

# Run models with different effects ---------------------------------------

# Effects for Jaccard out-degree, in degree popularity, out-in-degree 
# assortativity, and transitive triples type 1
m1_results <- estimate_saom(Jout, 100)
m2_results <- estimate_saom(inPop, 200)
m3_results <- estimate_saom(outInAss, 300)
m4_results <- estimate_saom(transTrip1, 400)

# Create a regression table -----------------------------------------------

htmlreg(l = list(m1_results, m2_results, m3_results, m4_results),
        omit.coef = "basic rate parameter dv",
        custom.coef.map = list("out-Jaccard similarity" = "Jaccard Similarity",
                               "indegree - popularity" = "In-degree Popularity",
                               "out-in degree^(1/2) assortativity" 
                                     = "Out-in-degree Assortativity",
                               "transitive triplets (1)" = "Transitive Closure",
                               "same aggression_var" = "Aggression Homophily",
                               "same subfaction_var" = "Subfaction Homophily",
                               "same militia_var" = "Militia Homophily",
                               "same role_var" = "Role Homophily",
                               "reciprocity" = "Reciprocity"),
         custom.model.names = c("Alliances", "Reputation", "Strong-vs-weak",
                                "Clustering"),
         caption = NULL,
         file = "figures/table1.html")

# Plot the results --------------------------------------------------------

plot_importance(m1_results, effect = "Out-degree\nJaccard similarity", 
                subtitle = "Alliance degradation")
ggsave("figures/figure3.tiff", width = 6.5, height = 4, dpi = 1000)

p7 <- plot_importance(m2_results, effect = "In-degree\npopularity", 
                      subtitle = "Preferential attachment") + 
   expand_limits(x = c(0, 0.6))
p8 <- plot_importance(m3_results, effect = "Out-in-degree\nassortativity", 
                      subtitle = "Aggressive vs non-aggressive attacks") + 
   expand_limits(x = c(0, 0.85))

p7 | p8

ggsave("figures/figure4.tiff", width = 8.5, height = 5.5, dpi = 1000)

plot_importance(m4_results, effect = "Transitive closure", 
                subtitle = "Clustering")
ggsave("figures/figure5.tiff", width = 6.5, height = 4, dpi = 1000)
